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Abstract 


An explicit-explicit subcycling procedure for the finite element analysis of structural 
dynamics is developed. This procedure has relaxed the usual constraint of requiring integer 
time step ratios for adjacent nodal groups. This allows for greater advantage to be taken of 
local stability criteria, and thus improves the efficiency of the explicit time integrator. 
Example problems are included to demonstrate the accuracy and stability of the method. 
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I. Introduction 


Explicit-explicit subcycling is oftentimes used in dynamic finite element problems 
in order to improve computational efficiency when the element size varies over the mesh. 
This is accomplished by separating the elements or nodes into groups and assigning a 
different time step to each group. The time step of each group depends only on the 
frequencies of the elements in that group. This results in great computational savings since 
less stiff elements can be assigned to the same group. In this way a time step is used for 
each group which is much closer to the critical time steps of the elements in the group. 
Mixed time integration was first introduced by Belytschko and Mullen[l] where they 
presented an "implicit-explicit" method for the second order equations which arise from the 
semidiscretization of the continuum equation. This was a nodal partition where both 
implicit and explicit methods were used in the time integration process. The stability of this 
type of procedure was studied in [2,3]- Hughes and Liu[4,5] later introduced an "implicit- 
explicit" method which divided the mesh into groups of elements rather than nodes, known 
as an element partition. Both methods were found to substantially improve computational 
efficiency. An "explicit-explicit" subcycling procedure was then introduced by Belytschko 
and Mullen[6]. In this procedure each group of nodes was integrated explicitly, but 
different time steps were used for each group. 

Mixed time integration procedures for first order semidiscretizations have also been 
proposed. Belytschko et al.[7] present sufficient conditions for stability of their element 
partition by using the element eigenvalue inequality directly on the evolution equations. 

The stability analysis of nodal partitions is complicated by the fact that their amplification 
matrices are unsymmetric. Smolinski et al.[8] avoid solving any unsymmetric systems in 
their nodal partition, and are able to develop stability conditions for semidiscretizations of 
the diffusion equation. 
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A major drawback to all these previous methods is that integer time step ratios are 
required. That is, if one group with a time step of At shares nodes or elements with a 
second group, then the only possibilities for the time step of the second group are [(n)At] 
or [At / n] where n is a positive integer. This introduces two difficulties when dividing the 
mesh into groups. First, it can severely restrict the division of the mesh into more than just 
a few groups, and therefore decrease the efficiency of the mixed time integration procedure 
being used. Second, since the time step is restricted to just a few possibilities, a time step 
may have to be used for a group that is much smaller than the stable time steps of the 
elements in the group. 

Recently, Mizukami[9] has relaxed this restriction for first order ordinary 
differential equations which result from semidiscretization of the diffusion equation. He 
presents an explicit-explicit nodal partition that allows for much greater flexibility when 
assigning time steps to the different nodal groups. He then gives a stability proof of the 
method, however, this proof only applies to a particular class of elements. Donea and 
Laval[10] have recently extended this proof of stability to arbitrary elements. 

The purpose of this paper is to present an explicit-explicit subcycling procedure 'for 
structural dynamics which does not require integer time step ratios for adjacent groups. It 
will be shown that this procedure requires very little additional storage and appears to be 
stable for a wide variety of example problems. 
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II. Nomenclature and definitions 


n 


At 


M 


At i 
Atf d 

At critical 


At e 
^ MLS 
1 M 

te 

C MLT 

d 


v 


a 

d e 


v e 


CT e 

a e 

e e 

e e 

f 

f ext 
f int 
ffnt 


= time step iteration counter 
= master time step 

= nodal time step 

= nodal time step for the previous step 
= critical time step of an element 
- time step used to update an element 
= master level time step 
= master clock 
= nodal clock 

= element clock 
= master level time 

= global vector of displacements 
= global vector of velocities 
= global vector of accelerations 
= element level vector of displacements 
= element level vector of velocities 
= Cauchy stress tensor for an element 

= stress rate tensor for an element 
= strain tensor for an element 

= strain rate tensor for an element 
= global vector of nodal forces 
= global vector of external nodal forces 
= global vector of internal nodal forces 
= element level vector of internal nodal forces 
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Wext 

= external work done on the system 

Wfcin 

= kinetic energy of the system 

Wi» t 

= internal energy of the system 

AWint 

= element level increment in internal energy 

B 

= gradient operator 

L e 

= element assembly operator 

M 

= lumped mass matrix 

mi 

= diagonal component of the mass matrix for node I 

Q e 

= domain of the element 

NUMEL 

= total number of elements 

NUMEL 

= total number of nodes 

MODOj) 

= remainder of the division of i by j 
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HI. Governing equations 


The equations of motion for structural dynamics are given by 

Ma + f = 0 (3.1) 

where, 

M = mass matrix 

a = global vector of accelerations 

f = global vector of nodal forces. 

The nodal forces are comprised of external and internal parts. 

f = (3-2) 

where, 

f ext = global vector of external nodal forces 

fi nt = global vector of internal nodal forces. 

The internal nodal forces for a single element are calculated from the element stresses and 
can be written 

ffnt = J n e BT<j e dQ (3.3) 


6 



where. 


f fnt = element internal force vector 

Q e = domain of the element 

B = gradient operator 

G e = Cauchy stress tensor for an element 

These internal forces are then assembled to the global vector 

fiat = fint+CLTffnt (3.4) 

where, 

L e = element assembly operator. 

For nonlinear analysis a strain rate formulation will be used to calculate the stress of 
an element. 


a e 

= C£ e 

(3.5) 

a e 

= <y e + At e 

(3.6) 


where, 

• ^ 

£ = strain rate tensor for an element such as the velocity strain or rate- 

of-defoimation 

• _ 

<7 = frame-invarient stress rate tensor for an element 
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C = constitutive matrix 

At e = time step used for the stress update of an element. 

The strain rate tensor is given by the strain displacement relation of the element as follows: 


e e = Bv e 


(3.7) 


where v e is the element velocity vector. 

In this paper, central difference time integration with variable time steps will be 
used to to advance the nodal velocities and displacements in time. The equations are given 
by the following: 


v I.n+ 1 /i ~ v I.n - 1 /2 + \ ( I + I Id ) a Ln 

d I.n+l = d I.n + I w l.n+ l t2 


(3.8) 

(3.9) 


where, 


a I.n 

= acceleration of node I at time (n)At 

v I.n + 1 /2 

= velocity of node I at time (n+V 2 )At 

d I.n+l 

= displacement of node I at time (n+l)At 

At i 

= time step of node I 

At° ld 

= time step of node I for the previous step. 


Since the time step for a node could change for each step, it is necessary to average 
the time steps as shown in equation (3.8) for the velocity update. Note that equation (3.6) 
has been written so that different time steps may be used for different elements. This 
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modification of the usual equations has also been made on equations (3.8) and (3.9) so that 
different time steps can be used for each node. The use of different time steps for different 
nodes and for different elements is common to all subcycling procedures. For the 
procedure to be presented here, however, the element update time steps of equation (3.6) 
are not necessarily the same as the nodal update time steps of equations (3.8) and (3.9). 
This added feature is unique to this new procedure and allows for the use of non integer 
time step ratios of adjacent nodes as will be explained in the following sections. 


9 



IV. Subcycling Procedure 


The basic idea underlying this subcycling procedure is the use of a master clock 
which is incremented by the smallest time step which occurs within the mesh; a master 
clock was also used by Belytschko and Liu[l 1]. This smallest time step will be referred to 
as the master time step. While the master clock is updated by the master time step, the 
nodal clocks are incremented by the time step of each node. Whenever a nodal clock time 
is equal to the master clock the velocity and displacement of the node are updated with 
equations (3.8) and (3.9) using the time step of the node. The time step of each node must 
be an integer multiple of the master time step and cannot be greater than the smallest critical 
time step of the elements surrounding the node. Although the nodal time steps must be 
integer multiples of the master time step, the attractive feature of this procedure is that it is 
not necessary for the time steps of two adjacent nodes to form an integer ratio. 

The times at which each element is updated are also governed by the nodal clocks of 
the nodes of the element Each element is updated every time one of the nodes of the 
element will undergo a velocity and displacement update. For the purpose of this element 
update, it is necessary to use an element clock for each element. This element clock is set 
equal to the master clock each time the element is updated. In this way a record is kept of 
when each element was last updated. The time step for the incremental update of the 
element stresses, equation (3.7), is then given by the difference between the master clock 
and the element clock. It is important to realize that this element update time step does not 
necessarily coincide with the time step of any of the element's nodes. 

As can be seen, the structure of this procedure is very similar to the subcycling 
algorithm described by Belytschko and Liu[l 1] and by Belytschko[12]. The major 
difference is that in the previous algorithm the times at which an element was updated was 
governed by element clocks while in this new procedure they are governed by the nodal 
clocks of the nodes of the element Figure 1 gives the nodal influence diagram for the 
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partition described in references [11] and [12]. In this figure an arrow represents the 
influence of one node on the velocity and displacement calculation of another node, a solid 
dot represents a time at which the velocity and displacement of a node are updated, and an 
open dot represents a time at which the displacement of a node would have to be 
interpolated between two known values. In order for a partition to work correctly the 
element stresses must be calculated whenever a node of the element is updated. In figure 1 
it can be seen that the stresses of element 1 are updated every 6At, while the stresses of 
elements 2 and 3 are updated every 2At as is shown by the bold horizontal lines. Because 
the time steps of these two nodes form an integer ratio, element stresses are calculated 
whenever they are needed for nodal updates. However, if this procedure is used in the 
case of non integer time step ratios as shown in figure 2, then the appropriate element 
stresses would not be available to update each node's acceleration, velocity, and 
displacement. For example, to update the nodal values of node 2 at time t = (n+3) At, it is 
necessary to have the element stresses at this time for all elements that are adjacent to node 
2. From figure 2 it is seen that the stress values of element 2 are known at times (n) At, 
(n+2) At, and (n+4) At, but not at the required time, (n+3) At If nodal clocks are used 
instead, as shown in figure 3, the values of the element stresses are calculated whenever 
they are needed for a nodal update so that even in the case of non integer time step ratios it 
is possible to correctly update the nodal values. 
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V. Implementation of the Subcycling Procedure 


A. Nodal Grouping 

The first step in the nodal grouping procedure is to assign a time step to 
each node. This is accomplished as follows: 

1. Assign a very large time step to each node, At j. 

2. Loop on the elements, (JE). 

Calculate the critical time step of the element, Ate nucal (JE). 

Loop on the nodes of this element, (I). 

At ! = Minimum(At l , At | ritical (JE)). ‘ (5.1) 

In this way each node is given the time step of the smallest element time step 
of the elements connected to the node. Next, it is necessary to determine a master 
time step so that all nodal time steps are an integer multiple of the master time step. 
This is most easily accomplished by letting the master time step be given by the 
smallest nodal time step. 

At M = Minimum ( At z ). (5.2) 

The remaining nodal time steps are then reduced by the equation 

Atj = At [Integer(At j / At )]. (5.3) 
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By equation (5.3) it is obvious that the master time step must be chosen so 
that it is smaller than or equal to all nodal time steps. It can also be seen from 
equation (5.3) that all nodal time steps are integer multiples of the master time step. 
Although equation (5.2) was used to determine the master time step, it should be 
noted that the master time step can be set to any value as long as the following two 
restrictions are satisfied. 


1. All nodal time steps are integer multiples of the master time step. 

2. A nodal time step is not greater than the smallest critical time step of the 
elements surrounding the node. 

Next, the master level time step must be calculated. This identifies the time 
increment that occurs between steps when all nodes are updated and is defined, 
therefore, as the least common multiple of all nodal time steps. 

At mls = Least common multiple ( At j ). (5.4) 


Once this value is calculated the master level time can be calculated by 
adding the master level time step to the current master time. 


1 MLT “ 1 M + MLT- 


(5.5) 


This master level time identifies the next time that all nodal and element 
clocks will be the same value, that is t j^. This means that all nodes and element 

stresses will have been updated to this master time level. It is useful to have all 
nodal values and element stresses updated to the same time level at some point in 
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the analysis so that the energy balance of the system may be checked-Once the 
master level time is determined the nodal grouping is complete. 

To illustrate the nodal grouping procedure, an example is given in figure 4, 
along with the element time steps, nodal time steps, master time step, master level 
time step, and master level time. 

B. Nodal Force Calculations 

After the master time step, nodal time steps, and the master level time have 
all been assigned as described in the previous section, the time stepping procedure 
can begin. The first step in the procedure is to calculate the nodal forces. The 
external nodal forces at each node are calculated as usual, however the calculation 
of internal nodal forces for each element must be modified. Before an element is 
updated the clocks of all the nodes of the element are checked by the condition 

Is (t x = t M ) for node I (5.6) 


where, 

tj = nodal clock 

t M = master clock. 

If condition (5.6) is true for any node of a given element, then the nodal 
clock is equal to the master clock and the node will be updated during this step. 
This nodal update will require the correct element stresses at this time for any 
element connected to this node. Therefore, this element must be updated during 
this step. This can be seen more easily with the help of an example. If it is 
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assumed that condition (5.6) is true for node I of an element, JE, then the velocity 
and displacement of node I will be updated by equations (3.8) and (3.9) which 
depend on the nodal acceleration, aj. The calculation of aj requires current internal 

nodal force values at node I, so that any element connected to node I must undergo 
a stress update at this time, so that the element stress coincides with the current 
time, t M . Therefore, if condition (5.6) is true for any node of an element, the 

stress of the element must be updated. For the purpose of this stress update an 
element time step must be calculated. This element time step is given by the 
equation 

At e = *M’*e (5.7) 

where, 

At e = element time step 

t e = element clock. 

This will give the time increment that has occurred since the last stress 
update of this element. Note that this element time step is not necessarily the same 
for an element each time it is updated. This can be seen by examining element 2 of 
figure 3. Once the element time step is calculated, the stress update of the element 

can take place. This update is given by equations (3.5) through (3.7) with the 
element update time step just calculated, At e . In this way the element is updated to 
the current master clock time, t M , using the appropriate time increment for each 

element. After the stresses have been updated the nodal forces are calculated in the 
usual manner by equations (3.3) and (3.4). The final step in the internal nodal 
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force calculation is to increment the clocks of all elements that have been updated 
with the equation 

te = ^ M' (5.8) 


C. Nodal Updates 


The next step in the time stepping procedure is to update the nodal variables. 
Again, central difference time integration is used, however, the procedure given by 
equations (3.1), (3.8) and (3.9) must be modified slightly. The calculation of the 
nodal acceleration at time (n)At is the same as before and is given by 


cfft - f ft> 


%n 


mi 


(5.9) 


for node I. To update the velocity of node I at time (n)At, the validity of the 
following condition is checked: 


Is (tj = t M ) for node I? (5.10) 

If condition (5.10) is true for node I, then equations (3.8) and (3.9) are 
used to update the velocity and displacement of node I. If condition (5.10) is false, 
then the nodal values are not changed. Once the new velocities and displacements 
have been calculated, the nodal clock is incremented by the nodal time step. 

t j — 1 1 *4" At j. (5 .11) 
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VI. Energy Calculations 


A major drawback to explicit time integration schemes is the fact that they are 
conditionally stable. When performing explicit calculations a limit must be placed on the 
time increment being used in order to prevent a loss of stability. In linear calculations this 
presents no problem since in this case the stability criteria remain the same for the duration 
of the analysis. However, in nonlinear calculations, the stability conditions of the mesh are 
constantly changing due to changes in problem geometry or material properties. Therefore, 
some means of checking the stability of the problem must be introduced. In linear 
problems, when an instability occurs it is characterized by exponential oscillatory growth in 
the nodal displacements. A periodic check on displacements can usually guarantee that no 
instability has occurred in the problem. However, in nonlinear situations where plastic 
deformations may dissipate large amounts of energy, the decreased modulus may enable 
the system to regain stability (known as an arrested instability, Belytschko[13]). In these 
situations simply checking for large unrealistic results may not always provide assurance 
that no instability has occurred, since an instability may be confined to a small region of the 
mesh and the problem may even regain stability at a later time. If this happens, the results 
obtained for the problem may seem quite reasonable, but the earlier instability in the 
problem will render the final results totally inaccurate. An alternative method for checking 
the stability of explicit finite element calculations is the use of an energy balance check. 
Hughes and Belytschko[14] gives this energy balance check in the form 

I W^+W^-Wr I < 5 HWlI (6.1) 


where, 


W„ xt = external work performed on the system up to time (n)At 
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W^ 11 = internal energy of the system at time (n)At 

Wn m = kinetic energy of the system at time (n)At 

8 = tolerance 

and IIWH is some measure of the total energy of the system. The use of any one of the 
terms on the left hand side of equation (6.1) alone would not provide a useful measure of 
this energy. For impulsively loaded problems when an initial velocity is prescribed, there 
would be no external work. For vibration problems the values of the kinetic energy and 
internal energy would oscillate and become very small at some point in the problem. 
Therefore, a measure of the total energy of the system should include all three terms. 

IIWH = IW^I+IW^I+IW^I. (6.2) 

The use of absolute value sign in the left hand side of equation (6.1) is not really 
necessary since an instability is usually characterized by a growth in kinetic and internal 
energies which are nonnegative except for roundoff errors; therefore, this left hand side 
would always be positive in the case of an instability. 

In time integration without subcycling equation (6.1) is checked at each time step to 
insure that the work done on the system and the energy of the system are nearly equal. 
Since a numerical instability in a problem always implies a loss of energy balance, it can be 
assumed that if equation (6.1) is true at every step, there is no such instability in the 
problem. For the case with subcycling, this energy balance cannot be checked at every 
step. This is because the internal energy of each element is not updated at every step. This 
can be seen in the equation for the increment of internal energy for an element. 

AW»‘-' = f^Jn* «J' n . m+ A<T=) T E'dn (6.3) 
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where. 


AWun ,e _ increment 0 f internal energy for an element from 
time (n-m) At M to (n) At M . 

In equation (6.3) the subscript m represents the number of master time steps that 
have passed since this element was last updated. After this increment of internal energy for 
an element is calculated, it is added to the total internal energy of the system. 

Wi" 1 = W^ + AW^ 1 ’ 6 . (6.4) 

Note that for the subcycling procedure, this element internal energy increment is not 
calculated for every element at every step. It is only calculated for an element when the 
stress of the element has been updated At other steps the internal energy in the element is 
not known. Therefore, the energy balance can only be checked at steps when the stress of 
every element is updated At these steps the contributions of all elements are added to the 
internal energy of the system and the correct value of the internal energy of the system is 
known. By definition these steps occur only at the master level times described earlier. 

The external work done on the system is calculated at every master time step and is 
given by the equation 

wr = W« 1 t +^-(vn- 1 /2)T(fext. f ext) > (6.5) 
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The kinetic energy can only be computed at half steps since the nodal velocities are 


only calculated at half steps. 

wSl h = 2 [( V n+1/2> T M v n+ l/ 2 ]. 


( 6 . 5 ) 


The kinetic energy at each step is then found by averaging the half step values. 

wf 1 = j(w„ k " l/j+ w n k 5 /2 ). (6.6) 
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VII. Additional Storage Requirements 

The addition of subcycling to an existing finite element code requires only four new 
additional storage requirements are as follows: 

a clock for each element 
a clock for each node 
nodal time step for each node 
nodal time step for each node for the previous step 
master level time step. 

The necessity of these additional storage requirements and a general overview of the 
procedure can be understood with the help of the flowchart given in figure 5. 


arrays of storage. The 


A tl 
At j ld 

^MLS 
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Vm. Numerical Examples 


To test the stability and accuracy of the proposed subcycling scheme, two example 
problems have been examined. Problem 1 consists of a bar subjected to a suddenly applied 
load at one end while the opposite end is fixed to a rigid wall. The bar is modelled with 32 
constant strain rod elements which have been divided into four groups. The problem 
dimensions and material properties are shown in figure 6, while the nodal grouping data is 
given in table 1. In case 1 the external force applied to the bar is small so that the problem 
re main s completely elastic, while in case 2 plastic deformation takes place due to the larger 
external force. The elastic problem is solved with the proposed subcycling procedure. 
Results of the elastic problem are given in figures 7 and 8. Figure 7 gives the time history 
of the stress at x=5.5 for both the subcycled and non subcycled cases. The exact solution 
is also included for reference. Figure 8 gives the energy error for both the subcycled and 
non subcycled cases. Figures 9 and 10 give similar results for the plastic problem. These 
results indicate that the subcycling procedure performs well for the proposed problem, 
since the stresses are similar with and without subcycling, and there is minimal energy 
balance loss for the subcycled case. 

The second example problem is a rotating ring subjected to a suddenly applied 
internal pressure. The problem dimensions and material properties are given in figure 1 1, 
and the finite element meshes are given in figures 12 and 14, Key[15]. Figure 13. gives 
the time history of the maximum stress of a single element for the mesh shown in figure 
12. Figure 15 compares the effective plastic strain of an element as a function of time with 
and without subcycling. Figure 16 gives the effective plastic strains of two different 
elements in the subcycled case, while the energy error is shown in figure 17. Figures 15 
through 17 all refer to the 100 element mesh of figure 14 with the nodal grouping data in 
table 2. From the figures it can be seen that the subcycled solution retains the symmetry of 
the problem and appears to be stable since there is no appreciable loss of energy balance. 
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VUI. Conclusions 


The subcycling, or mixed time integration, procedure described here has great 
flexibility in exploiting the increased stability of domains with larger of more flexible 
elements. Thus the difficulties associated with the presence of a few stiff elements are 
avoided. The procedure is easy to program, requiring on the order of 40 additional 
FORTRAN statements. Furthermore, it requires little additional storage: 

2 (NUMEL) + 2 (NUMNP) 


where, 


NUMEL = number of elements in the mesh 

NUMNP = number of nodes in the mesh. 

Stability proofs are not available for these subcycling procedures for second-order 
systems; stability proofs have been developed only for linear first-order systems, 
Belytschko, Smolinski, and Liu[7], Mizukami[9]. However, the energy stability checks 
clearly indicate that the procedure developed here is stable and we have not encountered any 
instabilities. 

Because of the great computational savings which can be achieved by this method, 
it is an attractive feature for any computer program used for engineering analysis by explicit 
time integration. 
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Table 1. Nodal grouping and time steps for example problem 1. 






Subcycling 

No subcycling 

Group 

Elements 

Nodes I 

ilsment tengtf 

i time step 

time step 

A 

1-10 

1-10 

1.0 

. lOAt 

At 

B 

11-15 

11-15 

0.4 

4At 

At 

C 

16-25 

16-26 

0.1 

At 

At 

D 

26-32 

27-33 

1.0 

10 At 

At 




Table 2. Nodal grouping and time steps for example problem 2. 
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Subcycling 

time step 
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time step 
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At 
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29-37 
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8At 

At 

3 

38-48 

73-94 

6At 

At 

4 

49-55 

95-108 

4At ' 

At 

5 

56-65 

109-128 

3At 

At 

6 

66-78 

129-154 

2At 

At 

7 

79-100 

155-200 

At 

At 
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Master to step 

At = 1.1 


Truncated nodal time steps 
Ati, Atg = 3At = 3.3 
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At 5 .At 6 .At 7 = At —1.1 
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Initialization 

1 . Read in data and set initial conditions, 

2 . Determine the values of At j, At M , and t MLS as described in section 
VILA. 

lime integration loop 

3 . Calculate the external nodal forces. 

4. Loop on the elements of the mesh (JE) 

If any node of element JE satisfies equation (5.5) 

THEN for element JE 

A. Calculate At M . with equation (5.6) 

B . Update the stress by equations (3.4) through (3.6). 

C. Calculate the internal nodal forces and assemble with 
equation (3.3). 

D. Calculate the increment in internal energy for this element and 
assemble to the internal energy of the system, equations (6.2) and (6.3). 

E. At e = At M 

5 . Use equation (1) to calculate the acceleration vector, ai. 

6. Loop on all the nodes of the mesh (I) 

A. If node I satisfies condition (5.9) 

THEN the velocity and displacement of node I are updated by equations 
(3.7) and (3.8). 

B. At? ld = At r 

C. If node I does not satisfy condition (5.9) 

THEN the velocity and displacement of node I are not updated. 

7 . Calculate the external work done of the system with equation (6.4). 

8 . Calculate the kinetic energy of the system with equations (6.5) through 

( 6 . 6 ) 

9 . If this is a master level time step, i.e., t M = t MLT 
THEN 

A. Check the energy balance with equation (6.1). 

B . If desired, recalculate the time steps, etc., as described in section VILA. 

10. Output the results. 

11. If the maximum time of the problem is not exceeded 
THEN 

A - 1 M = 1 M + At M' 

B. GOTO step 3. 

12. If the maximum time of the problem is exceeded 
THEN STOP. 
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